PHYS123 Group Homework Project¶

Wolf Mermelstein, Sara Elanchezhian, Mathew Gummere¶

December 2, 2023¶

Part G1¶

Calculate and plot the recursively defined logistics function $x_{n+1} = R x_n (1 - x_n)$.

We do this by creating a generator to yield logistics equation values, and then use it to find specific ranges of logistics equation values for given growth rates. For our samples we use 0.4 for $x_0$.

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from typing import Iterator

"""Generator that continuously yields logistics equation values at x_0 for some growth rate R."""
def compute_logistics_equation_values(R: float, x_0: float) -> Iterator[float]:
    x_prev = x_0
    while True:
        x_prev = R * x_prev * (1 - x_prev)
        yield x_prev
    
In [ ]:
"""Plot the logistic equation values for some growth rate R and initial value x_0."""
def plot_logistic_equation(iteration_counts, x_values, R, x_0):
    fig, ax = plt.subplots(figsize=(10, 6))

    # Create a line plot of the iteration counts vs. the logistic equation values.
    ax.plot(iteration_counts, x_values, marker='o', linestyle='-', color='b')
    ax.set_title(f'Logistic Equation (r={R}, x0={x_0})')
    ax.set_xlabel('Iteration')
    ax.set_ylabel('x')
    ax.grid(True)
    
"""
Compute the coordinates for the logistic equation values for some growth rate R and 
initial value x_0.
"""
def compute_logistics_equation_coordinates(R: float, x_0: float, N: int) -> np.ndarray:
    # Compute values using the generator that we previously defined.
    x_values = np.fromiter(compute_logistics_equation_values(R, x_0), float, count=N)
    # Create an array of iteration counts.
    iteration_counts = np.arange(N)
    return iteration_counts, x_values
    
# Plot the logistic equation values for some growth rate R and initial value x_0.
R_values = (1.5, 2, 2.8, 3.2, 3.6)
for R in R_values:
    x_coords, y_coords = compute_logistics_equation_coordinates(R, 0.4, 100)
    plot_logistic_equation(
        x_coords, 
        y_coords, 
        R=R,
        x_0=0.4
    )
    
plt.show()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Part G2¶

Finding values of $R_n$ where the logisics equation switches from

  • Period 1 to period 2, $R_1$
  • Period 2 to period 4, $R_2$
  • Period 4 to period 8, $R_3$
  • Period 8 to period 16, $R_4$
  • Period 16 to period 32, $R_5$
In [ ]:
"""
Generate bifurcation plot coordinates.

Arguments:
    equilibrium_value_generator: Generator that yields the equilibrium values given some R/G and x_0.
    value_count: Number of values to have for x axis.
    value_min: Initial value for x axis.
    value_max: Maximal value for x axis.
    iterations_to_start_at: Number of values to skip past for each x value.
    
Returns:
    The x and y coordinates.
"""
def compute_bifurcation_plot_coords(
        equilibrium_value_generator,
        value_count, 
        value_min, 
        value_max, 
        x_0=0.4,
        iterations_to_start_at=900, 
        x_per_value_count=100
    ):
    x_coords = np.zeros(value_count * x_per_value_count)
    y_coords = np.zeros(value_count * x_per_value_count)
    next_val = 0

    # For each R value in an evenly distributed range starting at <R_value_min> up to <R_value_max>
    # with R_value_count values overall.
    for x_type_val in np.linspace(value_min, value_max, value_count):
        # Create a generator to compute x values from
        generator = equilibrium_value_generator(x_type_val, x_0)
        
        # Skip past <iteration_to_start_at> number of computed x values
        for i in range(iterations_to_start_at):
            next(generator)
            
        # Generate <x_per_R_value_count> number of x values and add them to the plot at the current
        # R value (recall x axis is R)
        for i in range(x_per_value_count):
            x_coords[next_val] = x_type_val
            y_coords[next_val] = next(generator)
            next_val += 1

    return x_coords, y_coords

"""
Plot some x coordinates against y coordinates over a clipped region.

Args:
    x_coords: The x coords to plot.
    y_coords: The y coords to plot.
    title: The title of the plot.
    x_min: The minimum x cut off.
    x_max: The maximum x cut off.
    width: The width of the plot.
    height: The height of the plot.
    fontsize: The fontsize of the plot.
    dot_size: The size of the dots in the plot.
"""
def plot_diagram(
    x_coords, 
    y_coords, 
    title="",
    x_min=None, 
    x_max=None, 
    width=50, 
    height=30,
    fontsize=24,
    dot_size=0.05,
    gridline_count=30,
):
    # Create the plot
    fig, axes = plt.subplots(figsize=(width, height))

    # Clip the plot if we have a min or max provided
    if x_min or x_max:
        axes.set_xlim(x_min, x_max)

    # Plot the coordinates
    axes.scatter(x_coords, y_coords, marker='o', linestyle='-', color='b', s=dot_size)

    # Set the title and labels of the plot
    axes.set_title(f'{title} (from x={x_min} to x={x_max})', fontsize=fontsize)
    axes.set_xlabel('R')
    axes.set_ylabel('X_eq')
    
    # Specify the positions of our vertical gridlines by disabling the default x axis gridlines
    # and then manually adding in our own.
    axes.grid(axis='y', linestyle='--')
    gridlines = np.linspace(*axes.get_xlim(), gridline_count)
    axes.set_xticks(gridlines)
    for gridline in gridlines:
        axes.axvline(x=gridline, color='black', linestyle='-', linewidth=.3)
In [ ]:
# Compute the bifurcation plot coordinates
x_coords, y_coords = compute_bifurcation_plot_coords(
    compute_logistics_equation_values,
    20000, 
    2.5, 
    4,
)

"""An alias to plot a slice of the bifurcation diagram."""""
def plot_logistics_slice(x_min=None, x_max=None):
    plot_diagram(
        x_coords, 
        y_coords, 
        "Bifurcation Diagram", 
        width=300, 
        height=62,
        dot_size=.01,
        gridline_count=100,
        x_min=x_min,
        x_max=x_max,
    )
    
# Plot the bifurcation diagram and a slice of its tail so that later bifurcations are more
# clearly visible.
plot_logistics_slice()
plot_logistics_slice(3.5, 4)
    
plt.show()
No description has been provided for this image
No description has been provided for this image

Part G3¶

Repeat the exercise in HW G-2 but instead of using the logistics equation use the "laser" mapping $x_{n+1} = g x_n (1 - \tanh{(x_n)})$. Find the values of $G_n$ where the logistics equation switches from period 1 to period 2.

To achieve this we redefine a new generator, compute_laser_equation_values, to compute values, which we feed into the same x-y coordinate computing function, compute_bifurcation_plot_coords, that we defined earlier. We then plot various subsections to show period switches more clearly.

In [ ]:
from numpy import tanh

"""Generator that continuously yields laser equation values at x_0 for some growth rate R."""
def compute_laser_equation_values(G: float, x_0: float) -> Iterator[float]:
    x_prev = x_0
    while True:
        x_prev = G * x_prev * (1 - tanh(x_prev))
        yield x_prev
In [ ]:
# Compute the bifurcation plot coordinates for the laser equation
x_coords, y_coords = compute_bifurcation_plot_coords(
    compute_laser_equation_values,
    50000, 
    1, 
    10,
)

# Plot the entire bifurcation diagram
plot_diagram(
    x_coords, 
    y_coords, 
    "Laser Bifurcation Diagram, all periods", 
    width=100, 
    height=40,
    x_min=1,
    x_max=10,
    dot_size=.015
)

"""Plot a slice of the bifurcation diagram for the laser equation."""
def plot_laser_slice(x_min, x_max, title):
    plot_diagram(
        x_coords, 
        y_coords, 
        f"{title} Laser Bifurcation Diagram", 
        width=48, 
        height=26,
        x_min=x_min,
        x_max=x_max,
        dot_size=.04
    )
    
# Plot slices of the bifurcation diagram for the laser equation for different period switches
plot_laser_slice(5, 5.2, "Laser Bifurcation Plot, Period 0-1")
plot_laser_slice(8, 8.2, "Laser Bifurcation Plot, Period 1-2")
plot_laser_slice(9, 9.15, "Laser Bifurcation Plot, Period 2-3")
plot_laser_slice(9.2, 9.38, "Laser Bifurcation Plot, Period 3-4-5")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Part G5¶

Consider the logisics equation $x_{n+1} = R x_n (1 - x_n)$ with the growth parameter set to $R = 3.6$.

Define a sequence $x_{true}$ to be set of the $x_n$ calculated with $R = 3.6$ and $x_0 = 0.5$

Define a sequence $x_{pred}$ to be set of the $x_n$ calculated with $R = 3.6$ and $x_0 = 0.5 + \epsilon $.

The idea is that even a small error $\epsilon$ in the initial conditon will lead to predicted values that vary wildly from the "true" values.

In [ ]:
"""Add an error to an x_0 value."""
def with_error(x_0, error):
    return x_0 + error

"""Build an array of x values for a given error value."""
def compute_x_values(x_0, error=0):
    return np.fromiter(compute_logistics_equation_values(R, with_error(x_0, error)), float, count=N)

"""Superpose an additional plot onto the main plot axes."""
def plot_y_axis(axes, x_values, y_values, color, label, width=.5):
    axes.scatter(x_values, y_values, color=color, label=label, linewidth=width, s=.1)
    axes.set_ylabel(label, color=color)

# Create the plot
fig, axes = plt.subplots(figsize=(40, 16))
plt.title("Effect of x_0 errors on bifurcations")

# Create y axis for the various errors
errors = (0, 0.05, 0.005, 0.0005)
colors = ('black', 'blue', 'red', 'green')

# Plot the data
for error, color in zip(errors, colors):
    x_coords, y_coords = x_coords, y_coords = compute_bifurcation_plot_coords(
        compute_logistics_equation_values,
        3000, 
        2.5, 
        4,
        x_0=with_error(x_0, error),
    )
    label = f'Error = {error}'
    plot_y_axis(axes, x_coords, y_coords, color=color, label=label)
    
# Add a legend
axes.legend(loc='upper left', fontsize=28)
axes.grid(True)

# Show the plot
plt.show()
No description has been provided for this image
In [ ]: